Has the ITV changed at high elevation across two study watersheds in Colorado, and if so, how? 1) Determine the temperature bias correction for individual SNOTEL stations within the study sites. 2) Compute ITV after removing seasonal and long-term climate change-induced trends. 3) Examine the elevational, seasonal, and spatial distribution of ITV across the two study watersheds.
There are three bias correction methods to address the first objective: the National Oceanic and Atmospheric Administration (NOAA) 9th order polynomial correction (or “NOAA9”) recommended by the Snow Survey (Atwood et al., 2023), the 4th order polynomial correction based on the co-located new and old sensor sites in Idaho (referred to as the “Morrisey” method) (Ma et al., 2019), and the Oyler et al., (2015) dataset.
Start new document
NOAA6 is better and not as unrealistically complicated as NOAA9
Trend differences? Statistically different? Among the GHCN and adjacent non-GHCN station
Examine station metadata – how often and much (?) did they move
“write everything down, these are our options, what we want to do, what we don’t”
Accumulation (start of accumulation to peak), melt (peak to SAG), snow free (SAG to start of accumulation)
Fix the period, or average start of accumulation/melt/SAG?
Are these the same dates for all stations, or the same for each station
Yes, the only reason we didn’t originally try it was due to many SNOTEL minimums being -51°C, which was erroneous.
“Is it better than what Oyler did? Low elevation stations don’t represent higher elevation stations.” Need to discuss what we don’t do.
We need to at least discuss this.
As Tmax, Tmin, with code in hand, this should be easy
Working again with Crosho:
SNOTEL_426 <- snotel_download(site_id = 426, path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_426,"C:/Users/steen/OneDrive/Documents/R/Projects/zombie_code/data_raw/snotel_426.csv", row.names = FALSE) #write in the raw data
snotel_426 <- read.csv("C:/Users/steen/OneDrive/Documents/R/Projects/zombie_code/data_raw/snotel_426.csv", header = TRUE)
Morrisey 7/21/2005
#str(snotel_426) # check the date, usually a character.
snotel_426$Date <- as.Date(snotel_426$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_426_clean <- snotel_426 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_426_clean <- snotel_426_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_426_clean, aes(x = Date, y = temperature_mean)) +
geom_point(alpha =0.3) + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
# # For SF average annual temp
# sf <- snotel_426_clean %>%
# group_by(waterYear) %>%
# mutate(anntemp = mean(temperature_mean)) %>%
# distinct(anntemp)
#
# ggplot(sf, aes(x=waterYear, y= anntemp))+
# geom_point()+
# geom_smooth(method = "lm", se=TRUE)+
# theme_few()
Raw SNOTEL 426 average temperatures for water years 1986-2021
snotel_426_clean <- snotel_426_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -50) %>%
filter(temp_diff < 40)
ggplot(snotel_426_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature difference (°C)') +
xlab('Date')
Maximum minus minimum SNOTEL 426 temperatures for water years 1986-2021
# filtering for too few observations in a year
snotel_426_cull_count_days <- snotel_426_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_426_cull_count_days
## # A tibble: 3 × 2
## # Groups: waterYear [3]
## waterYear n
## <dbl> <int>
## 1 1988 281
## 2 1989 340
## 3 2009 340
snotel_426_clean_culled <- snotel_426_clean %>%
filter(waterYear != "1988" & waterYear != "1989" & waterYear != "2009")
# ggplot(snotel_426_clean_culled, aes(x = Date, y = temp_diff)) +
# geom_point() + #lwd = 2) +
# theme_few() +
# ylab('Daily temperature varience (°C)') +
# xlab('Date')
ggplot(snotel_426_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)
dygraph(temp_426_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
Cleaned SNOTEL 426 average temperatures for water years 1986-2021
Looking at minimum and maximum values
# filtering for too few observations in a year
min_temp_426_xts <- xts(snotel_426_clean_culled$temperature_min, order.by = snotel_426_clean_culled$Date)
dygraph(min_temp_426_xts) %>%
dyAxis("y", label = "Daily min temperature (°C)")
SNOTEL 426 minimum temperatures for water years 1986-2021
max_temp_426_xts <- xts(snotel_426_clean_culled$temperature_max, order.by = snotel_426_clean_culled$Date)
dygraph(max_temp_426_xts) %>%
dyAxis("y", label = "Daily max temperature (°C)")
snotel_426_clean_culled <- snotel_426_clean_culled %>%
mutate(temp_difference = (temperature_max-temperature_min))# %>%
#filter(temperature_mean > -50) %>%
#filter(temp_diff < 40)
SNOTEL 426 maximum temperatures for water years 1986-2021
difference_temp_426_xts <- xts(snotel_426_clean_culled$temp_difference, order.by = snotel_426_clean_culled$Date)
dygraph(difference_temp_426_xts) %>%
dyAxis("y", label = "Temperature difference (max-min) (°C)")
SNOTEL 426 temperature difference (max-min) for water years 1986-2021
# ggplot(snotel_426_clean_culled, aes(x = Date)) +
# geom_point(aes(y = temperature_min, color = "min"), size =0.7, alpha =0.2)+
# geom_point(aes(y = temperature_max, color = "max"), size =0.7, alpha =0.2)+
# geom_vline(xintercept = specific_date, linetype=4, size=0.7)+
# theme_few() +
# geom_smooth(aes(x=Date, y = temperature_min, color = "min"), method=lm, se=FALSE)+
# geom_smooth(aes(x=Date, y = temperature_max, color = "max"), method=lm, se=FALSE)+
# scale_colour_manual(name = "daily temps", values=c("min"="blue", "max" = "red"))+
# ylab('Temperature (°C)') +
# xlab('Date')
#Basic:
# ggplot(snotel_426_clean_culled, aes(x = Date)) +
# geom_point(aes(y = temperature_min, color = "min"), size =0.7, alpha =0.2)+
# geom_point(aes(y = temperature_max, color = "max"), size =0.7, alpha =0.2)+
# geom_vline(xintercept = specific_date, linetype=4, size=0.7)+
# theme_few() +
# geom_smooth(aes(x=Date, y = temperature_min, color = "min"), method=lm, se=FALSE)+
# geom_smooth(aes(x=Date, y = temperature_max, color = "max"), method=lm, se=FALSE)+
# scale_colour_manual(name = "daily temps", values=c("min"="blue", "max" = "red"))+
# ylab('Temperature (°C)') +
# xlab('Date')
before_subset <- filter(snotel_426_clean_culled, Date <"2005-07-21")
after_subset <- filter(snotel_426_clean_culled, Date >= "2005-07-21")
ggplot() +
geom_point(data = before_subset, aes(x = Date, y = temperature_min, color = "min"), size =0.7, alpha =0.2) +
geom_smooth(data = before_subset, aes(x = Date, y = temperature_min), method = "lm", se = FALSE, color = "blue") +
# Plot the second subset with another trendline
geom_point(data = after_subset, aes(x = Date, y = temperature_min, color = "min"), size =0.7, alpha =0.2) +
geom_smooth(data = after_subset, aes(x = Date, y = temperature_min), method = "lm", se = FALSE, color = "blue")+
geom_point(data = before_subset, aes(x = Date, y = temperature_max, color = "max"), size =0.7, alpha =0.2) +
geom_smooth(data = before_subset, aes(x = Date, y = temperature_max), method = "lm", se = FALSE, color = "red") +
# Plot the second subset with another trendline
geom_point(data = after_subset, aes(x = Date, y = temperature_max, color = "max"), size =0.7, alpha =0.2) +
geom_smooth(data = after_subset, aes(x = Date, y = temperature_max), method = "lm", se = FALSE, color = "red")+
geom_vline(xintercept = specific_date, linetype=4, size=0.7)+
theme_few() +
scale_colour_manual(name = "daily temps", values=c("min"="blue", "max" = "red"))+
ylab('Temperature (°C)') +
xlab('Date')
Daily SNOTEL 426 minimum and maximum temperatures for water years 1986-2021
Crosho 426 Morrisey 7/21/2005
Need four datasets:
# 2) NRCS pre-sensor change adjusted data (Morrisey) *"morrisey"*
snotel_426_adjusted <- snotel_426_clean_culled %>%
mutate(morrisey = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean)) %>%
mutate(noaa_conus = if_else(Date >= "2005-07-21", 610558.226380138*(((temperature_mean+65.929))/194.45)^9-2056177.65461394*(((temperature_mean+65.929))/194.45)^8+2937046.42906361*(((temperature_mean+65.929))/194.45)^7-2319657.12916417*(((temperature_mean+65.929))/194.45)^6+1111854.33825836*(((temperature_mean+65.929))/194.45)^5-337069.883250001*(((temperature_mean+65.929))/194.45)^4+66105.7015922199*(((temperature_mean+65.929))/194.45)^3- 8386.78320604513*(((temperature_mean+65.929))/194.45)^2+ 824.818021779729*(((temperature_mean+65.929))/194.45)-86.7321006757439, temperature_mean)) %>%
mutate(noaa_morrisey = 610558.226380138*(((morrisey+65.929))/194.45)^9-2056177.65461394*(((morrisey+65.929))/194.45)^8+2937046.42906361*(((morrisey+65.929))/194.45)^7-2319657.12916417*(((morrisey+65.929))/194.45)^6+1111854.33825836*(((morrisey+65.929))/194.45)^5-337069.883250001*(((morrisey+65.929))/194.45)^4+66105.7015922199*(((morrisey+65.929))/194.45)^3- 8386.78320604513*(((morrisey+65.929))/194.45)^2+ 824.818021779729*(((morrisey+65.929))/194.45)-86.7321006757439)
#Adding min-max
snotel_426_adjusted <- snotel_426_adjusted %>%
mutate(min_morrisey = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_min^(4))+(3.72*10^(-5))*(temperature_min^(3))-(2.16*10^(-3))*(temperature_min^(2))-(7.32*10^(-2))*(temperature_min)+1.37)+temperature_min, temperature_min)) %>%
mutate(noaa_conus = if_else(Date >= "2005-07-21", 610558.226380138*(((temperature_min+65.929))/194.45)^9-2056177.65461394*(((temperature_min+65.929))/194.45)^8+2937046.42906361*(((temperature_min+65.929))/194.45)^7-2319657.12916417*(((temperature_min+65.929))/194.45)^6+1111854.33825836*(((temperature_min+65.929))/194.45)^5-337069.883250001*(((temperature_min+65.929))/194.45)^4+66105.7015922199*(((temperature_min+65.929))/194.45)^3- 8386.78320604513*(((temperature_min+65.929))/194.45)^2+ 824.818021779729*(((temperature_min+65.929))/194.45)-86.7321006757439, temperature_min)) %>%
mutate(noaa_min_morrisey = 610558.226380138*(((min_morrisey+65.929))/194.45)^9-2056177.65461394*(((min_morrisey+65.929))/194.45)^8+2937046.42906361*(((min_morrisey+65.929))/194.45)^7-2319657.12916417*(((min_morrisey+65.929))/194.45)^6+1111854.33825836*(((min_morrisey+65.929))/194.45)^5-337069.883250001*(((min_morrisey+65.929))/194.45)^4+66105.7015922199*(((min_morrisey+65.929))/194.45)^3- 8386.78320604513*(((min_morrisey+65.929))/194.45)^2+ 824.818021779729*(((min_morrisey+65.929))/194.45)-86.7321006757439)
#max
snotel_426_adjusted <- snotel_426_adjusted %>%
mutate(max_morrisey = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_max^(4))+(3.72*10^(-5))*(temperature_max^(3))-(2.16*10^(-3))*(temperature_max^(2))-(7.32*10^(-2))*(temperature_max)+1.37)+temperature_max, temperature_max)) %>%
mutate(noaa_conus = if_else(Date >= "2005-07-21", 610558.226380138*(((temperature_max+65.929))/194.45)^9-2056177.65461394*(((temperature_max+65.929))/194.45)^8+2937046.42906361*(((temperature_max+65.929))/194.45)^7-2319657.12916417*(((temperature_max+65.929))/194.45)^6+1111854.33825836*(((temperature_max+65.929))/194.45)^5-337069.883250001*(((temperature_max+65.929))/194.45)^4+66105.7015922199*(((temperature_max+65.929))/194.45)^3- 8386.78320604513*(((temperature_max+65.929))/194.45)^2+ 824.818021779729*(((temperature_max+65.929))/194.45)-86.7321006757439, temperature_max)) %>%
mutate(noaa_max_morrisey = 610558.226380138*(((max_morrisey+65.929))/194.45)^9-2056177.65461394*(((max_morrisey+65.929))/194.45)^8+2937046.42906361*(((max_morrisey+65.929))/194.45)^7-2319657.12916417*(((max_morrisey+65.929))/194.45)^6+1111854.33825836*(((max_morrisey+65.929))/194.45)^5-337069.883250001*(((max_morrisey+65.929))/194.45)^4+66105.7015922199*(((max_morrisey+65.929))/194.45)^3- 8386.78320604513*(((max_morrisey+65.929))/194.45)^2+ 824.818021779729*(((max_morrisey+65.929))/194.45)-86.7321006757439)
# ggplot(snotel_426_adjusted, aes(x = Date)) +
# geom_point(aes(y = temperature_min, color = "min"), size =0.7, alpha =0.2)+
# geom_point(aes(y = temperature_max, color = "max"), size =0.7, alpha =0.2)+
# geom_vline(xintercept = specific_date, linetype=4, size=0.7)+
# theme_few() +
# geom_smooth(aes(x=Date, y = temperature_min, color = "min"), method=lm, se=FALSE)+
# geom_smooth(aes(x=Date, y = temperature_max, color = "max"), method=lm, se=FALSE)+
# scale_colour_manual(name = "daily temps", values=c("min"="blue", "max" = "red"))+
# ylab('Temperature (°C)') +
# xlab('Date')
adj_before_subset <- filter(snotel_426_adjusted, Date <"2005-07-21")
adj_after_subset <- filter(snotel_426_adjusted, Date >= "2005-07-21")
ggplot() +
geom_point(data = adj_before_subset, aes(x = Date, y = noaa_min_morrisey, color = "min"), size =0.7, alpha =0.2) +
geom_smooth(data = adj_before_subset, aes(x = Date, y = noaa_min_morrisey), method = "lm", se = FALSE, color = "blue") +
# Plot the second subset with another trendline
geom_point(data = adj_after_subset, aes(x = Date, y = noaa_min_morrisey, color = "min"), size =0.7, alpha =0.2) +
geom_smooth(data = adj_after_subset, aes(x = Date, y = noaa_min_morrisey), method = "lm", se = FALSE, color = "blue")+
geom_point(data = adj_before_subset, aes(x = Date, y = noaa_max_morrisey, color = "max"), size =0.7, alpha =0.2) +
geom_smooth(data = adj_before_subset, aes(x = Date, y = noaa_max_morrisey), method = "lm", se = FALSE, color = "red") +
# Plot the second subset with another trendline
geom_point(data = adj_after_subset, aes(x = Date, y = noaa_max_morrisey, color = "max"), size =0.7, alpha =0.2) +
geom_smooth(data = adj_after_subset, aes(x = Date, y = noaa_max_morrisey), method = "lm", se = FALSE, color = "red")+
geom_vline(xintercept = specific_date, linetype=4, size=0.7)+
theme_few() +
scale_colour_manual(name = "daily temps", values=c("min"="blue", "max" = "red"))+
ylab('Temperature (°C)') +
xlab('Date')
NOAA9 over Morrisey corrected SNOTEL 426 minimum and maximum temperatures for water years 1986-2021
426 temperature_mean mean SD and median SD NON-CORRECTED.
#average water year temperature
yearly_wy_aver_426 <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_426 <- yearly_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(temperature_mean))
#average mean temperature by day for the period of record:
daily_wy_aver_426 <- daily_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_426$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_426 <-daily_wy_aver_426 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_426$date_temp <- (daily_wy_aver2_426$date_temp) #reduce the sig figs
# ggplot(daily_wy_aver2_426, aes(x = waterDay, y = date_temp))+
# geom_line(size= 0.7) +
# theme_few() +
# ylab('Average Daily temperature (°C)') +
# xlab('Day of water year')
# mean SD
standard_dev_426 <- daily_wy_aver_426 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426 <- standard_dev_426 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426 <- standard_dev_all_426 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 3.824933 |
| 1990 | 3.918921 |
| 1991 | 4.205106 |
| 1992 | 3.875164 |
| 1993 | 3.884094 |
| 1994 | 3.913071 |
| 1995 | 4.051010 |
| 1996 | 4.040079 |
| 1997 | 4.113160 |
| 1998 | 3.669836 |
| 1999 | 4.089543 |
| 2000 | 3.843753 |
| 2001 | 3.838346 |
| 2002 | 4.144204 |
| 2003 | 3.914148 |
| 2004 | 4.121643 |
| 2005 | 3.659180 |
| 2006 | 4.098573 |
| 2007 | 3.931999 |
| 2008 | 3.898354 |
| 2010 | 3.813336 |
| 2011 | 3.841523 |
| 2012 | 3.620427 |
| 2013 | 4.080408 |
| 2014 | 3.926416 |
| 2015 | 3.934119 |
| 2016 | 3.644881 |
| 2017 | 4.367449 |
| 2018 | 3.767432 |
| 2019 | 3.755422 |
| 2020 | 4.085168 |
| 2021 | 3.826197 |
| 2022 | 4.103787 |
# adding median SD
standard_dev_all_426_med <- standard_dev_426 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426_med <- standard_dev_all_426_med %>%
group_by(waterYear) %>%
mutate(resid_median = median(residual)) %>%
mutate(sd_1_med = residual-resid_median) %>%
mutate(sd_2_med = (((sum((sd_1_med)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2_med, .keep_all = TRUE) %>%
select(waterYear, sd_2_med)
standard_dev_all_426_med %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2_med |
|---|---|
| 1987 | 3.832605 |
| 1990 | 3.918993 |
| 1991 | 4.217058 |
| 1992 | 3.908350 |
| 1993 | 3.936908 |
| 1994 | 3.948249 |
| 1995 | 4.055615 |
| 1996 | 4.075571 |
| 1997 | 4.123535 |
| 1998 | 3.676261 |
| 1999 | 4.095045 |
| 2000 | 3.846280 |
| 2001 | 3.865732 |
| 2002 | 4.190361 |
| 2003 | 3.921409 |
| 2004 | 4.125068 |
| 2005 | 3.666076 |
| 2006 | 4.106820 |
| 2007 | 3.964453 |
| 2008 | 3.916529 |
| 2010 | 3.817120 |
| 2011 | 3.841673 |
| 2012 | 3.620881 |
| 2013 | 4.133183 |
| 2014 | 3.944373 |
| 2015 | 3.934234 |
| 2016 | 3.644881 |
| 2017 | 4.372454 |
| 2018 | 3.768991 |
| 2019 | 3.777040 |
| 2020 | 4.205946 |
| 2021 | 3.835441 |
| 2022 | 4.104938 |
standard_dev_all_426 <- standard_dev_all_426 %>%
left_join(standard_dev_all_426_med, by= 'waterYear')
ggplot(standard_dev_all_426, aes(x = waterYear))+
geom_point(aes(y = sd_2, color = 'average'), size =2, alpha =0.4)+
geom_smooth(aes(y = sd_2, color = 'average'), method = lm, se= FALSE)+
geom_point(aes(y = sd_2_med, color = 'median'), size =2, alpha =0.4)+
geom_smooth(aes(y = sd_2_med, color = 'median'), method = lm, se= FALSE)+
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Mean and median standard deviation of SNOTEL 426 average temperatures for water years 1986-2021
426 temperature_mean mean SD and median SD noaa_morrisey CORRECTED.
#using the clean culled df:
#average water year temperature
yearly_wy_aver_426_noaa_morrisey <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_noaa_morrisey = mean(noaa_morrisey))
#Average temperature by day for all water years:
daily_wy_aver_426_noaa_morrisey <- yearly_wy_aver_426_noaa_morrisey %>%
group_by(daymonth) %>%
mutate(aver_day_temp_noaa_morrisey = mean(noaa_morrisey))
#average mean temperature by day for the period of record:
daily_wy_aver_426_noaa_morrisey <- daily_wy_aver_426_noaa_morrisey %>%
group_by(daymonth) %>%
mutate(all_ave_temp_noaa_morrisey = mean(daily_wy_aver_426_noaa_morrisey$aver_day_temp_noaa_morrisey))
# try to show all years as means.
daily_wy_aver2_426_noaa_morrisey <-daily_wy_aver_426_noaa_morrisey %>%
group_by(waterDay) %>%
mutate(date_temp_noaa_morrisey = mean(noaa_morrisey))
daily_wy_aver2_426_noaa_morrisey$date_temp_noaa_morrisey <- signif(daily_wy_aver2_426_noaa_morrisey$date_temp_noaa_morrisey,3) #reduce the sig figs
# ggplot(daily_wy_aver2_426_noaa_morrisey, aes(x = waterDay, y = date_temp_noaa_morrisey))+
# geom_line(size= 0.7) +
# theme_few() +
# ylab('Average Daily temperature (°C)') +
# xlab('Day of water year')
standard_dev_426_noaa_morrisey <- daily_wy_aver_426_noaa_morrisey %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_noaa_morrisey-aver_ann_temp_noaa_morrisey)+noaa_morrisey-aver_day_temp_noaa_morrisey) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426_noaa_morrisey <- standard_dev_426_noaa_morrisey %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426_noaa_morrisey <- standard_dev_all_426_noaa_morrisey %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426_noaa_morrisey %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 3.830349 |
| 1990 | 3.917975 |
| 1991 | 4.175217 |
| 1992 | 3.900360 |
| 1993 | 3.886360 |
| 1994 | 3.836302 |
| 1995 | 4.095875 |
| 1996 | 4.061514 |
| 1997 | 4.130496 |
| 1998 | 3.652421 |
| 1999 | 4.162786 |
| 2000 | 3.879900 |
| 2001 | 3.794383 |
| 2002 | 4.051192 |
| 2003 | 3.887463 |
| 2004 | 4.156371 |
| 2005 | 3.730215 |
| 2006 | 4.347486 |
| 2007 | 4.199739 |
| 2008 | 4.164094 |
| 2010 | 4.104876 |
| 2011 | 4.027864 |
| 2012 | 3.896582 |
| 2013 | 4.321041 |
| 2014 | 4.122428 |
| 2015 | 4.117491 |
| 2016 | 3.877980 |
| 2017 | 4.611047 |
| 2018 | 4.017755 |
| 2019 | 4.019420 |
| 2020 | 4.387524 |
| 2021 | 4.145641 |
| 2022 | 4.343523 |
# ggplot(standard_dev_all_426_noaa_morrisey, aes(x = waterYear, y = sd_2))+
# geom_line(size= 0.7) +
# theme_few() +
# geom_smooth(method = "lm", se=FALSE) +
# ylab('SD') +
# xlab('Water year')
# adding median SD
standard_dev_all_426_med_noaa_morrisey <- standard_dev_426_noaa_morrisey %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426_med_noaa_morrisey <- standard_dev_all_426_med_noaa_morrisey %>%
group_by(waterYear) %>%
mutate(resid_median = median(residual)) %>%
mutate(sd_1_med = residual-resid_median) %>%
mutate(sd_2_med = (((sum((sd_1_med)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2_med, .keep_all = TRUE) %>%
select(waterYear, sd_2_med)
standard_dev_all_426_med_noaa_morrisey %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2_med |
|---|---|
| 1987 | 3.831038 |
| 1990 | 3.918001 |
| 1991 | 4.187841 |
| 1992 | 3.912562 |
| 1993 | 3.918909 |
| 1994 | 3.849165 |
| 1995 | 4.096066 |
| 1996 | 4.073020 |
| 1997 | 4.135222 |
| 1998 | 3.656442 |
| 1999 | 4.175637 |
| 2000 | 3.886817 |
| 2001 | 3.821890 |
| 2002 | 4.088239 |
| 2003 | 3.893692 |
| 2004 | 4.156372 |
| 2005 | 3.730336 |
| 2006 | 4.355936 |
| 2007 | 4.248006 |
| 2008 | 4.202192 |
| 2010 | 4.108711 |
| 2011 | 4.036565 |
| 2012 | 3.901007 |
| 2013 | 4.391629 |
| 2014 | 4.149276 |
| 2015 | 4.117681 |
| 2016 | 3.884976 |
| 2017 | 4.623524 |
| 2018 | 4.023378 |
| 2019 | 4.058160 |
| 2020 | 4.533625 |
| 2021 | 4.174225 |
| 2022 | 4.350527 |
standard_dev_all_426_noaa_morrisey <- standard_dev_all_426_noaa_morrisey %>%
left_join(standard_dev_all_426_med_noaa_morrisey, by= 'waterYear')
ggplot(standard_dev_all_426_noaa_morrisey, aes(x = waterYear))+
geom_point(aes(y = sd_2, color = 'average'), size =2, alpha =0.4)+
geom_smooth(aes(y = sd_2, color = 'average'), method = lm, se= FALSE)+
geom_point(aes(y = sd_2_med, color = 'median'), size =2, alpha =0.4)+
geom_smooth(aes(y = sd_2_med, color = 'median'), method = lm, se= FALSE)+
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
NOAA9 over Morrisey corrected mean and median standard deviation of SNOTEL 426 average temperatures for water years 1986-2021